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CHAPTER I 
INTRODUCTION 

In many non-linear and/or time-varying problems, solutions are 
often sought for the "best" or optimal way to control or guide a system 
subject to some criterion. Typically, this criterion consists of an 
objective fr "tional (pay-off function) that is to be extremized 
(maximum or minimum), the desired terminal conditions, and the possible 
constraints on the state and/or control. After the problem has been 
mathematically formulated, there are various conceptual techniques 
available to obtain the conditions required to numerically solve the 
optimization problem. Most o f these conceptual techniques may be 
considered in one of two main categories, direct or indirect. Direct 
methods are dependent on the direct evaluation of an objective 
functional in order to arrive at control changes, which in turn, result 
in an improvement in the objective functional. Indirect methods, on 
the other hand, use variational analysis to arrive at a set of 
conditions that h e optimal solution must satisfy without actually 
evaluating the objective functional. Thus, the problem becomes one of 
finding solutions which satisfy the conditions of a resulting two-point 
boundary value (TP3V) problem. 

The most attractive of the direct methods are the gradient methods 
developed by Kelley [1]* simultaneously with those of Bryson and 


♦Numbers in brackets refer to the references in the Bibliography at the 
end of this thesis. 
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and Denham [2], These methods have the advantage that convergence is 
not contingent upon good initial guesses and that improvement in the 
criterion to which it is subject is guaranteed by each step. The 
methods have the great disadvantage that convergence often 
deteriorates in the 'icinity of the optimum, since the gradient goes 
to zero. However, Lasdon, Mitter, and Warren [3] have used second 
derivative information to improve the convergence near the optimum. 
Numerous other extensions have been made to the basic gradient 
methods, and these are lucidly discussed by Sage [4] and Wilde and 
Beightler [5], 

Among the more usual indirect approaches are the calculus of 
variations as explained among others by Bryson and Ho [6] and Bliss [7]; 
the dynamic programming of Dreyfus and Bellman [8,9]; and Pontryagin's 
Principle [10]. A disadvantage of the majority of these methods is 
that a solution is often contingent upon good initial guesses, 
particularly when complex problems are involved. The;e methods, 
however, do possess the great advantage that whenever a solution is 
obtained, its associated control is considered to be optimal because 
of the necessary conditions which must be satisfied to obtain a solution. 
The calculus of variations concept, in particular, offers a powerful 
tool for solving optimization problems because of the conditions 
resulting from this classical theory. 

This thesis presents a solution to a complex lifting reentry 
three-degree-of-freedom problem by using the concept of the calculus of 
variations to generate a set of necessary conditions. Then solving 
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numerically the resulting TPBV problem by an improved existing 
technique and a new "direct-indirect" technique, Both of these 
techniques are designed to decrease the strong dependency of the 
solution upon good initial guesses. Since such mathematical treatment 
leads to a TPBV problem whose solution requires usage of large memory 
high-speed computers, the solution of practical problems has been 
rather limited until the last decade. 

The atmospheric entry of a lifting body is a formidable physical 
problem. When optimal solutions are sought to entry problems, several 
simplifying dynamic and kinematic assumptions must be made if results 
are to be obtained with a reasonable amount of computational effort. 

This study addresses itself to the motion of the center of gravity of 
the vehicle as it passes through the atmosphere of a non-rotating 
spherical earth. The only external forces considered are the gravita- 
tional and the dissipative aerodynamic forces; vehicle control is 
effected through lift vector modulation. An approximate expression 
due to Detra, Kemp, and Riddell [11] is used to calculate the stagna- 
tion convective heat rate, neglecting other heat transfer processes. 

The three-dimensional trajectory formulation is tailored for ease of 
numerical solution and proper representation of the problem at hand. 
Previous entry studies conducted by the author and vacuum trajectory 
studies by Tapley, Szebehely, and Lewallen [12] and Lewallen, Schwausch, 
and Tapley [13] show that effort in formulating differential equations 
for numerical solutions often pays dividends by way of increased accu- 
racy and speed. For this reason, the equations of motion describing 
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the entry of a lifting body are in a spherical coordinate system and 
are non-dimensionalized. In most dissipative force problems such as 
1 this one, the medium is an important factor in the resulting motion. 

| With this in mind, the atmosphere was not modeled in the conventional 

I manner shown in the literature (such as in Chapman [14], Loh [ 1 5], and 

Klafin and Barnhard [16]) p = p 0 exp (-h/h Q ), but was modeled by a 

J '{+] 

logarithmic form such as lnp = a Q +2 a^h . In this manner, a very 

good model of the atmosphere is obtained that represents a continuous 
atmosphere between any desired two altitudes. 

In many problems in the field of flight mechanics, it is 
important to restrict the total heat input to a vehicle during entry 
into a planetary atmosphere along with the total aerodynamic load. 

The optimization problem may now be stated as follows: G'Wen a set 

of differential equations which describe the motion of a lifting 
vehicle entering through the earth's atmosphere, find the combination 
of state and control variables "'hat will minimize the integral on time 
of the sum of the aerodynamic force and the heat flux (pay-off function), 
and at the same time, satisfy specified terminal conditions which 1 ill 
determine the terminal trajectory time. As posed, this problem falls 
into a class of control optimization problems due to Bolza [17]. The 
entry problem considered here does not have state and/or control con- 
straints along the trajectory. Problems with such constraints are 
fully discussed by Bryson and Ho [6], Hestenes [18], and Lastman and 
Tapley [!"]. Proceeding with the calculus of variations on the posed 
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problem results in a TPBV problem with the appearance of a set of 
co-states (Lagrange multipliers) which have little physical 
significance and often are a discouraging factor in optimization 
studies. This is due to the fact that an initial guess of the co-states 
must be made to start the numerical integration of the differential 
equations associated with the TPBV problem. More often than not, the 
result is an initial solution which is totally unacceptable from both 
the numerical viewpoint and the physical viewpoint. Thus, it should 
be obvious that the optimal solution of a complex problem is dependent 
upon good initial guesses and that it follows that a major effort of 
this study is to attempt to desensitize the dependence of an optimal 
solution on the initial guesses. 

Having formulated the TPBV problem, it now remains to choose a 
method to solve it. Specifically, the problem at hand is to 
systematically meet the boundary conditions while estimating the 
initial co-states. Numerical techniques to accomplish this hava been 
considered by Hestenes [20] as early as 1949 and have been improved 
steadily since then. Some of tne popular methods include the method 
of adjoint functions of Jarovics and McIntyre [21] which uses the 
equations adjoint to the linearized co-state equations; the 
quasilinearization methods formulated through a concept as presented 
by Kalaba [22] along with the modified quasi 1 inearizati on method by 
Lewallen [23]; the perturbation methods based on the work of Breakwell 
[24] and Breakwell et al [25]; and, more recently, the method of 
perturbation functions (MPF) of Lewallen [23]. Based on the work of 
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Lewallen [23], which compares the convergence characteristics of the 
MPF, the modified quasilinearization method, and the method of 
adjoint functions, among others on an orbital transfer problem, the 
author chose the MPF because of its reported superior performance. 

The MPF method was modified to improve the iteration philosophy 
for correcting the terminal dissatisfaction (terminal boundary errors). 
The iteration philosophy of Lewallen [23] was based on a norm reduction 
method only. More recently, Lastman and Tapley [19] have compared the 
MPF method using a norm reduction method and a magnitude limitation 
method separately with limited success. Also recently Doiron [47] has 
improved the MPF method by essentially employing a correction factor 
based on the norm and the desired number of iterations. In this thesis, 
the author uses a combination norm reduction and magnitude limitation 
method defined as an improved method of perturbation functions (IMPF). 

A new method which uses one of three function minimization techniques 
on a terminal dissatisfaction function during the early iterations and 
, then switches to the IMPF is also exercised on the entry problem. This 

l 

I 

i method is referred to as the modified method of perturbation functions 

(MMPF). The three function minimization techniques used are the pattern 

i 

! search method of Hooke and Jeeves [26], Davi don's [27] conjugate 

| gradient search method of Fletcher and Powell [28], and an accelerated 

random search technique [29]. In order to test the ability of the IMPF 
and the three MMPF techniques to converge when started with different 
I guesses, the entry optimization problem was started with six different 

sets of guesses. 
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There have been several optimization techniques applied to the 
problem of reentry trajectories with partial success. Bryson et al [30] 
used the gradient method to solve a planar entry problem for minimum 
entry heat, and Breakwell et al [25] solved a planar entry problem 
for maximum terminal speed at a given altitude. The Pontryagin 
maximum principle was used by Leondes and Niemann [31] on a planar 
entry problem. Lastman and Tapley [19] considered a two-degree-of- 
freedom entry using an objective functional similar to the one used 
in this study and employing the MPF technique. Tapley and Williamson 
[32], along with Colunga [33], considered several Apollo entry problems 
using MPF with an objective functional similar to the one used here. 

No where in the literature has the author been able to find an entry 
problem using the same features as the one presented in this thesis; 
one that could be used directly as a basis of comparison for the IMPF 
and MMPF strategies. Therefore, an MPF strategy similar to Lewallen's 
[23] is used on several different initial guesses to serve as a 
comparison basis. 

In the following chapter, the entry problem is completely defined 
from a general statement of the problem through the formulation of the 
pay-off function. Chapter III presents the formulation of an optimiza- 
tion problem using the calculus of variations to arrive at a set of 
conditions which guarantee optimality. In addition, the resulting 
TPBV problem is defined. In Chapter IV a solution of a TPBV problem 
by the MPF is outlined that includes the differential correction of 
norm reduction and magnitude limitation which give rise to the IMPF. 
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The MMPF using each of the three function minimization schemes on a 
terminal dissatisfaction function is presented in Chapter V. Chapters 
III through V address an optimization problem in general. Chapter VI 
specifically focuses on the entry problem stated in Chapter II and 
presents its boundary conditions along with a discussion and analysis 
of the attempt to desensitize the effect of initial guesses on con- 
vergence. Finally, a short summary and the conclusions of this study 
appear in the last chapter. 
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CHAPTER II 
PROBLEM DEFINITION 

2.1 Statement of the Entry Problem 

The statement of the problem for the present study is to 
find the path which a spacecraft, using roll modulation for control 
with given aerodynamic characteristics and entering through the 
earth's atmosphere, must follow in going from one point in state 
space to another with a certain objective functional. The objective 
functional requires that the total integral of the sum of the 
stagnation convective 'neat flux and t h s total aerodynamic forces 
experienced by the vehicle during its trajectory between the two 
points in space be minimized, at least in a local sense. The 
objective function, as stated above, results in a design requirement 
which favors the spacecraft thermal protection system weight and 
exposes the crew to short-duration, high-level aerodynamic force pulses. 

2.2 Basic Assumptions 

The aerodynamic properties, mainly the lift and drag 
coefficients (C^, Cp), are assumed to remain constant for a given 
angle of attack in the hypersonic flight regime [34], For purposes 
of this study, the angle of attack was assumed constant; consequently, 
trajectory control is available through roll angle modulation only. 

The gravitational field is an inverse square, gravitational force 
field, and the earth and its atmosphere are non-rotating with respect 
to inertial space. 
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In most dissipative force problems, such as the one being 

presented, the dissipative medium is an important factor in the ensuing 

motion. With this in mind, the atmosphere [35] was not modeled in the 

classical manner shown in the literature, i.e. p = p Q exp (-h/h Q ), but 

j -i+1 

was modeled by a logarithmic form such as lnp = a„ + 2 a. h . The 

0 i=l 1 

advantage of choosing the latter functional form over the classical 
form is that a continuous atmosphere may be defined over any altitude 
range to any desired degree of accuracy rather than a piecewise con- 
tinuation. Thus, a very accurate model of an exponential atmosphere 
was generated by obtaining the "a" coefficients from a least squares 

j_ L. 

fit t; the order of an actual atmosphere [36], For this study, 
j=6 was fou nH to be sufficiently accurate. Taking the exponential 
of the logarithmic form mentioned earlier, the exponential density, 


p = S exp 


6 

2 

i=0 


(e-jh 1 ) 


results. 


2.3 Governing Differential Equations of Motion 

The mathematical model used in deriving the differential 
equations of motion assumes, in addition to the assumptions in 
Section 2.2, a point mass with three degrees-of-freedom whose motion 
is referenced to an inertial X, Y, Z coordinate system and is expressed 
in a spherical coordinate system (r, y, e) as shown in Figure 2-1. 
Control of the spacecraft is effected by rolling the stability axis 
about the velocity vector (Figure 2-2), thus using the vehicle's lift 
to perform out-of-plane maneuvers (roll modulation). The two 


4 


! 


I 





Figure 2-2. Vehicle Lift, Drag, and Control Angle Definition 
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aerodynamic lift force components resulting from a roll maneuver u, 
are defined as 

L y * p V 2 AC l cos u (2-la) 

L h = \ p V 2 AC l sin u (2-lb) 

where -it u it and the aerodynamic drag force as 

D = ipV 2 AC D (2-lb) 

Then using the inverse square gravitational field 

r 2 

9 = 9 0 {■—) (2-ld) 

and vehicle mass 

tn = jf- (2-le) 

y o 

it may be show' [37,38] that the six first-order differential equations 
describing the motion are 

»- -Jp r * % ( r- )Z si " r (2 - 2a > 


r = -V sin r (2-2b) 

V 2 C l A g 0 Sgp ♦ (^£) 2 g 0 S2SJ1 ( 2 - 2c ) 

1 ’ J» V 2 C l A jjy ? ^ s u r g 0 - ^ cot 6 cos r sin A (2-2d) 


0 

Y 


V cos r cos a 

r 

V cos r sin A 
r sin 9 


(2-2e) 

(2-2f) 
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The above set of equations (2-2) form an n set of non-linear 
differential equations whose solution in general is obtained by a 
numerical integration process. Numerical solutions are also quite 
dependent on the range of values which the dependent as well as the 
independent variables assume. Studies presented in the literature 
suggest that a normalization scheme of some sort aids in the solution 
of this type of equations. Therefore, the procedure used to non- 
dimensional ize the equations of motion (Equation 2-2) involved deriva- 
tives as well as dependent variables. The transformation applied in 
going from dimensional to non-dimensional derivatives was 


d( ) _ r e 

-zr-T c 


41 

ir 


where 


V = /g r 
c ^o e 

In addition, the following length, velocity, density, and energy 
constants were used to non-dimensionalize the dependent variables: 

Length - Radius of the earth, r g 20.925738 x 10® ft 

2 

Gravity - Sea level gravity, g Q 32.174 ft/sec 

Density - Sea level density, p Q 0.2378 x 10 ^ slug/ft^ 

Energy - Mechanical to thermal energy 

conversion factor, J 778.0 ft-lb/Btu 

Thus, the set of n equations (2-2) normalized according to the above 
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scheme with a scalar control u are as follows: 

_ 1 ' :,2 C D A , sin r 

V = - J p V -r- + — 

c W r 


(2-3a) 


r = - V sin r 

r = - \ cos r - i p v C. cos i* 

r £ L W vr 


(2-3b) 

(2-3c) 


a' = - l cos r - \ o V C, sin u A + 

r * u W Mr 


(2-3d) 


V cos r co c a 


(2-3e) 


/ = v c ? s r sin A (2-3f) 

r sin e 

2.4 Formulation of the Objective Functional Argument 

As mentioned in Section 2.1, the objective functional is the 
minimum integral of the sum of the convective stagnation heat flux and 
the aerodynamic forces evaluated between the initial and final 
trajectory points. The argument of this objective functional is then 
the heat flux and the aerodynamic forces. The convective stagnation 
heat flux expression commonly used in ths literature to des ribe the 
heat input to a spacecraft is that due to Detra, Kemp, and Riddel [11]. 
Referenced to a one-foot radius sphere, the heat flux can be 
expressed as 


A * 17 fion ( - — -.3.15 / p \i 

h i/,oju i 26>000 ; ' 0 . 002378 ' 
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and the total aerodynamic load as 

n 2 L v 2 l h 2 

G - (g) + (g~) + (g-) (2-5) 

Normalizing Equations (2-4) and (2-5) as specified in 
Section 2.3, the following is obtained: 

Q = C 1 V 3 * 15 p 0,5 (2-6) 

G = i l (C 2 + C 2 ) 1/2 ^V 2 (2-7) 

c U 

= Q + pG (2-8) 

where p is a weighing parameter, and ij; is the argument of the 
objective functional. 

The problem as defined in this chapter may be classified as 
an optimization problem whose necessary conditions and their 
applications are derived and shown respectively ir. Chapter III. 

2.5 Specified Initial and Terminal Conditions 

As pointed out in Section 2.1, the problem under study is 
to find the path a spacecraft must follow to satisfy some objective 
functional as described in Section 2.4 while going from some specified 
initial state to a specified final state. The specified initial state 
is defined as the initial conditions of the problem and the final state 
<<c the terminal conditions. For the problem under study, the initial 
boundary conditions specified are n+1 ; That is, all of the independent 
variables n as well as the dependent variable time, which is zero, are 
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specified. The terminal conditions specified are all of the indepen- 
dent variables except the heading (azimuth) A; thus, at the terminal 
time there are n-1 specified terminal conditions. 


1 


CHAPTER III 


FORMULATION OF AN OPTIMIZATION PROBLEM 

3.1 Necessary Conditions for Optimal Trajectories 

Trajectory optimization problems classically require the 
satisfaction of various conditions. These evolve from consideration 
of the following problem: Given a set of non-linear differential 

equations* determine the history of its variables so that some 
objective functional is optimized while satisfying specified initial 
and terminal constraints. The objective functional is, in general, 
some combination of the variables describing the problem. Variables 
in optimization problems are classically divided into state, co-state, 
independent (normally time), and control variables. 

The non-linear set of differential equations which describe 
the trajectory are 

X = f (x, u, t) (3-1) 

where, for the problem under study, x is an n-vector of state 
variables; u is an m-vector of control variables; and t is the 
independent variable, time. The initially specified initial 
conditions are 

L <x 0 , t 0 ) = 0 (3-2) 

where L is a p-vector, and the initially specified terminal conditions 
with M as a q-vector are 


M (x f , t f ) = 0 


(3-3) 
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The objective function^, the quantity to be extremized, is a 
scalar and, in general, is formulated as a problem oi Bolza [17] 

}f 

P ■ $ (X f , t f ) +/ ip (x, u, t) dt (3-4) 

t o 

where 4» = «(> (x^, t^) allows the Introduction of functions whose 
terminal values must be extremized. The non-linear differential 
equation (3-1) may be adjoined to P with the co-state X: 

P = 4> + / L + X T (f - x)l dt (3-5) 

*0 

For convenience, the scalar term H is defined as follows: 

H = \p + X^ f (3-6) 

where H is commonly referred to as the variational Hamiltonian 
function. The term H takes its name from the well known "Hamiltonian 
function" of classical mechanics [39] because of its functional 
similarity. Substituting Equation (3-6) into (3-5) results * : n the 
following form of P: 

tf 

P = <j > + / (H - X T x) dt (3-7) 

t o 

Equations (3-1) through (3-3) and (3-7) describe in general 
many optimization problems for deterministic, non-linear, time- 
dependent systems. 
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Integrating the second term under the integral of the above 
expression (3-7) by parts, then 


P = <j> - A x 


tf T 

t + «{ (H + A X 

0 0 


) dt 


(3-8) 


The first necessary condition classically required 
throughout the literature [6] for an optimal trajectory by using 
calculus of variations techniques is that the lirst variation of P 
vanish. Citron [40] shows that in th° calculus of variations theory, 
the first variation is identical to the total differential in 
optimization of functions. The conditions are discussed by Tapley 
and Lewallen [41] and a derivation is included for completeness. 
Taking the first variation of P in the form of (3-8), the result is 
given by 


dP = d<t> 


f 


- (qa T X + A T dx) 


(H + A x) dt 


(3-9) 


Taking the indicated differentials in the above expression and using 
the Liebnitz rule on the last term. Equation (3-9) becomes 


dP 


= (<f>^ dx + <j>£ dt) 1^ - (dA^ x + A^ dx) 

f u 

t f t 

+ (H + a T x) dt t + 1 f [(H l + a T f x ) fix + hJ 6 u 
+ 6a + A T fix + fiA T xj dt (3-10) 
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Performing the indicated integration on the last term in 
(3-10) and defining the total differential of \ to the first order as 

T T T T 

dx = 6X + x dtp where the variation in x, 6X , takes place during 
a fixed time period, the necessary condition that dP vanish now 


becomes 


dP = - X T ) dx t + (<t> t + H) dt t + (x T dx - Hdt) J t 

+ { f [(hJ + x T ) 6x + hJ 6u + <SX T (f - x)] dt = 0 


The condition that dP vanish implies that each term in 
(3-11) goes to zero provided that all of the variations are 
independent. The resulting necessary conditions may be classified 
as initial transversal ity conditions, terminal transversal ity 
conditions, and conditions to be satisfied for all time t between t 


and tp. 


The initial transversal ity condition 


X 1 dx - H dt , =0 (3-1 

% 

will be identically satisfied if the initial state and time are 
specified. However, if the initial state and time are not specified 
and dx and dt are independent, then (3-12) implies that x T and H must 
vanish at t Q . On the other hand, if dx and dt are not independent, 
use must be made of (3-2) to satisfy (3-12). 
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The remaining two transversal ity conditions which must be 
satisfied at the terminal time t f are 


and 


+ H) dt 



(3-13) 


(*J " xT) dx |t f = 0 


(3-14) 


The above two equations imply that if dt and dx are not specified, 
their respective coefficients, (<t> t + H) and (4J - x T ), must be zero. 

The only remaining necessary conditions derived from dP = 0 
(3-11) come from within the integral term in that expression and must 
be satisfied for all t Q < t < t^. The first condition within the 
integral of (3-11) is 


hJ + X T = 0 (3-15) 

which provides '■•he co-state differential equations and is the 
classical Euler Lagrange equation. The next condition is the classical 
optimality condition 

H u = 0 (3-16) 

The third and last condition resulting from the integrand of (3-11) 
is merely that the original differential equations of motion be 
satisfied, that is 

x - f = 0 (3-17) 
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A trajectory satisfying the necessary conditions expressed 
by Equations (3-13) through (3-17) defines a stationary trajectory, 
these conditions are not sufficient to guarantee an optimal trajectory. 
A fourth condition, and one that is su^icient to insure an optimum, 
is known as the Legendre-Clebsch condition in the calculus of 
i variations. This condition involves the Weierstrass E-Function as 

explained by Gelfand and Fomin [42] as well as Bryson and Ho [6] and 
must be equal to, or greater than, zero for a minimum. The 

■w 

Legendre-Clebsch condition states that, for t < t < t^, H u(j must be 
non-singular positive definite 


H uu > 0 (3-18) 

3.2 Reduction of an Optimization Problem to a Two-Point Boundary 
Value Problem 

A trajectory optimization problem may now be reduced to a 

TPBV problem by considering all of the conditions for optimality 

derived in Section 3.1. The conditions that must be satisfied for 

all t < t < t- are as follows: 
o f 

x = f = h[ (3-19) 

where x is an n vector first-order, non-linear, differential equation 
of motion; 

A = - H l (3-20) 

where A is also an n vector first order non-linear differential 
equation of co-states and is called the Euler-Lagrange equation; 

H u = 0 (3-21) 
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which is the classical optimality condition. Since u is an m vector, 
(3-21) defines m algebraic equations. These m algebraic equations 
may be solved for m control variables, u, as a function of the state, 
co-state, and independent variable time and thus u may be eliminated 
from Equations (3-19) and (3-20) if desired; 


V,* 0 


(3-22) 


this inequality states that H uu must be non-singular positive 
definite. It is commonly used in trajectory optimization problems 
to resolve sign ambiguities arising out of the condition for u 
from the optimality condition H u = 0. 

Equations (3-19) and (3-20), the state and co-state 
equations, may be re-written in a more compact form as 


= F (Z,u,t) (3-23) 


where Z is a 2n vector composed of state and co-state variables. 

The initial boundary conditions which must be satisfied at 



• 

X 


(x, A ,u ,t) 




T 


A 


1 

t 

m 

X 

X 

* 

* 


t = t are 
o 


9 (x Q , t Q ) = 0 (3-24) 

which is an n vector made up of the initially specified initial 
conditions L from Equation (3-2) and the initial transversal ity 
conditions from Equation (3-12). As pointed out in Section 3.1, if 
dx and dt in Equation (3-12) are not independent, which is generally 
the case, then dL = 0 is used to eliminate p variations from the n+1 
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variations of Equation (3-12), leaving n+l-p independent variations 
whose coefficients are equated to zero to yield n+l-p relations. Thus, 
when these n+l-p relations are combined with the p relations arising 
out of L, there results a set of n+1 relations, g = 0 and t Q . 

h (Z f , t f ) = 0 (3-25) 

is an n+1 terminal boundary condition vector which is made up of the 
initially specified terminal conditions M from Equation (3-3) and the 
terminal transversal ity conditions from Equations (3-13) and (3-14). 

An argument parallel to that used above to insure that dx and dt are 
independent on the initial boundary also applies here on the terminal 
boundary. Thus, n+1 independent coefficients are assured. 

Equations (3-24) and (3-25) completely define the 2n+2 
conditions necessary to solve a TPBV problem. In many problems, the 
initial time t Q is usually specified as zero; therefore, under this 
assumption, we have defined n initial conditions, g, and n+1 terminal 
conditions, h. In summary, the TPBV problem to be solved is defined 
by a system of 2n non-linear, first-order, differential equations (3-23) 
with n initial boundary conditions, g, and n+1 terminal boundary 
conditions, h. The solution of this problem results in a minimum of 
the objective functional (3-4) with satisfaction of all the constraints 
defined by (3-2) and (3-3). The general solution of these problems is 
so complex that closed-form solutions are only realized under a 
multitude of simplifying assumptions, if at all; thus, solutions must 
be obtained numerically. In fact even numerical solutions are complex. 
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The following two chapters will discuss this problem and use 
perturbation methods to obtain solutions. Chapter IV will deal with 
an improved version of the popular method of perturbation functions 
( IMPF ) and Chapter V will address a new approach, a modified method 
of perturbation functions (MMPF). Chapter IV also outlines an 
algorithm employing the improved version of MPF which, with minor 
modification, is also applicable with the MMPF outlined in Chapter V. 


i 



CHAPTER IV 


SOLUTION OF A TPBV PROBLEM BY AN IMPROVED 
METHOD OF PERTURBATION FUNCTIONS 


4.1 Introduction 

Equations (3-21. - ; through (3-25) are the transformation of 
the entry problem as stated in Chapter II into an optimal control 
problem. This resulting optimal control problem is a TPBV problem 
whose numerical solution may be obtained in general only through 
some sophisticated iteration algorithm. The need for an iteration 
algorithm arises from the fact that there are 2n non-linear, first- 
order differential equations (3-23) which must be numerically 
integrated either forward or backwards in time; however, there are 
only either n or nfl conditions known for the forward and backwards 
integration, respectively. Thus, it is necessary to guess either n 
or n+1 conditions to allow integration of the differential equations. 

An improved method of perturbation functions (IMPF) is used to update 
the assumed conditions and to continue iteratively the updating 
procedure until the terminal boundary conditions h go to zero. 

Figure 4-1 shows the general outline of such an algorithm. 

4.2 Method of Perturbation Functions 

As the name itself implies, the method consists of perturbing 
the Z function. This function, Equation (3-23), with the control u 
eliminated by using Equations (3-21) and (3-22) is perturbed linearly 
about the k^ trajectory at the end of each iteration. This pertur- 
bation then provides an estimate of the value of the unknown variables 
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needed to begin the Integration process at the next iteration and, at 
the same time, drive the dissatisfaction, h, to zero. In general terms, 
this process is Illustrated In Figure 4-1. Practical use of this 
method has shown that it possesses second nrder convergence properties 
[41] when it is near a solution 

If the initial boundary conditions satisfy g - C, 

■ 

Equation (3-24), then Z - F (Z,t). Equation (3-23), may be integrated 
forward in time starting at t * 0 to seme assumed t^. Generally, 
h U on the initial iteration sc that Equation (3-25) is not 
satisfied, mainly because the n assumed co-states along with the 
final time estimate, t f , are not the optimal values. The terminal 
dissatisfaction function h is then used tc correct the assumed values 
by the procedure described below. 

Consider a perturbed trajectory terminating on a perturbed 
set of terminal conditions. If only the linear terms of a Trylor Series 
expansion about some reference trajectory are retained, then 

* 5 k < It l k - V <•-’> 

and if 6Z - Z£ + -j - then the above results can be expressed as 

6Z = |£ | k <5Z (4-2) 

which is a set of 2n linear pertur^ tion equation with being a 

x L. 

2n x 2n matrix evaluated on the k n or reference trajectory. Since, 
on the reference trajectory h f 0, a perturbation equation is also 
obtained for the terminal constraint h which is evaluated at the 
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terminal time tf. The differential constraint dh (Z^,t^) is 

dh = Iz I tf dZ f + ll lt f dt f (4 ‘ 3) 

and represents a change in the terminal error due to changes in 
and t f . If dZ f = 5Z f + Z f dtf is used to allow terminal Z variations 
due to tf, then Equation (4-3) becor v js 

dh = <5Zf + h dt^ (4-4) 

The desired linear algebraic correction equation is derived 
under the conditions that there exists a fundamental matrix n, as 
defined in the literature, so that the terminal variations in Z (6Zf) 
can be related to the initial variations in Z (sZ Q ) according to 

<5Z f = tt (t f ,t Q ) 6Z q (4-5) 

The above expression is the well-known solution to a differential 

equation of the form cf Equation (4-2). Equation (4-5) implies that 

6 Z f is some linear function of the initial variation 6Z„, where it is 
t o 

a 2n x 2n matrix. Thus, Equation (4-4) can now be written in a more 
useful form as 

dh = fj | t * 6Z q + h dt f (4-6) 

by making use of Equation (4-5). As stated earlier, n of the 2n 
elements in SZ Q are known; therefore, no corrections are required for 
these states. Thus, only the right half of the 2n x 2n ir matrix in 
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Equation (4-5) needs to be known, and it may be rewritten as 


5Z f = *JX 5i O 


(4-7) 


where ir^ x is a 2n x n matrix. Using this reduced form. Equation (4-6) 
in partitioned matrix form now becomes 


dh 




(4-8a) 


or, upon solving for the correction vector and assuming that T is 
non-singular. 


--S] = T' 1 dh (4— 8b ) 

dt f J 


There remains two variables to be determined in the above 
expression, namely, and dh. Since the estimate of dh will be 
covered in the following section, the determination of can be 
made. Evaluation methods for ir or tt,, are well covered in the 

OA 

literature [43] but will be summarized here for the sake of 
completeness. Consider Equation (4-5), which is a solution to 
Equation (4-2), and relates the final perturbations to the 
initial perturbations 6Z Q . Assume that a 2n x 2n identity matrix I 
is used to express the initial perturbations where each column of I 
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is used as in the following manner: 


fl . 0 0 

‘ 1,0 
0 1 

' 0 , 


io 

10 


I = 


1 ~ 

L° 


0 

0 

0 


(4-9) 


where i = 1 ,2n refers to each column of I. Then each successive 

XL 

application of 6Z results in a <$Z f which is an i tn column of the 
transition (fundamental) matrix ir. Therefore, if Equation (4-2) is 
integrated 2n times with 6Z Q as an initial vector from t Q = 0 to t^, 

the transition matrix tt would be available. Furthermore, for the 
case at hand where only is needed, the integration of Equation (4-2) 
is required only n times, and the initial vector used is Sl Q where 
a = n+1 , 2n. 

4.3 Determination of the Differential Correction 

The usual procedure found in the literature (Lewallen [23] 
is an example) for the evaluation of dh is to let 


dh = - n h (4-10) 

with the restriction that the sca^r correction factor 0< n < 1 and 
the values that n then assumes are chosen as some step function of 
llh^. Equation (4-8b) may now be rewritten by substituting for dh with 
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Equation (4-10) as 


* - nT _1 h (4-11) 

For problems whose initial solutions result in large ||h|[, n is 
normally chosen close to zero; as the norm of the dissatisfaction ||h|| 
begins to vanish r approaches one. The successful solution to a 
problem often depends to a large degree on the logic used for 
selecting n. 

Thus, the logic used to develop a process for selecting n 
should concentrate on the fact that an attempt is being made to solve 
a highly complex non-linear problem by at best a quasi -second order 
technique. This means that particular attention should be directed 
toward the following: 

a) The case where n = 1 which indicates 100 percent of the 

[6X I 

correction vector — - is requested. Large corrections 

[dt f J 

should in general be suspect and not allowed. 

b) The maximum (100 percent) correction should be based 
on the expected range of values of the state. The 
comparison can be made by comparing the norms of the 
correction vector and the range of values of the state. 

c) The norm of the terminal dissatisfaction should also be 
monitored. Steps to reduce large percentages of h 
should not be allowed basically because elimination of 
large errors in a non-linear system which is based on 
linear estimates has little guarantee for success. 
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The general philosophy of the correction factor using the 
logic stated in (a), (b), and (c) above is to avoid overshooting the 
envelope of convergence; improve the area of convergence; and, at 
the same time, speed up the convergence process. It will later be 
shown in Chapter VI that this philosophy does, in fact, improve the 
envelope of convergence as compared to results in the literature. 

The functional relationship for the correction factor is 


■f 

( 1 / 


(4-12) 


(ilhl!) if ddilS Hell 

11/ MH if ||d ij > |je II 

where f (||hi|) is defined as a step function, lldll is the norm of 100 
percent of the correction vector (Equation (4-11) evaluated with 
n = 1), and Hell is a constant representing the norm of the expected 
range of the states and time of the particular problem whose solution 
is being sought. In the event that the present norm is greater than 
the past, n is reduced by a factor of two. 

4.4 Iteration Scheme for the Improved Method of Perturbation Functions 


Now that the initial correction vector equation 


SA. 


i_dt f . 


= -nT" 1 h 


has been completely defined, an iteration scheme for solving the TPBV 

problem may be stated. The scheme requires that an initial vector of 

A and an estimate of the terminal time, t^, be made available so that 
* 

the 2n equations, Z = F(Z,t), may be integrated forward in time to 

t = t f . Simultaneously with this integration, the 2n perturbation 
* dF 

equations 6Z = 6Z are integrated to the same t^ n times using the 
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starting vector defined by the right half of Equation (4-9). At the 
terminal time, the terminal dissatisfaction h as well as its norm, is 
evaluated. Subsequently, the algebraic equations which determine the 
initial corrections in x and t^ to be applied at the next iteration 
(Equation 4-11) are solved. This Iteration procedure is continued 
until some small predetermined value of Ithll is reached, when it is 
assumed that the problem is solved or has converged. However, it is 
not always possible to achieve convergence under ordinary circumstances. 
If convergence is not achieved within a reasonable number of iterations 
and divergence is diagnosed, it is best to start again with a different 
estimate on the unknowns. 

The described iteration scheme for the IMPF may be sequen- 
tially applied in the following manner: 

a) Choose initial estimates for x q , t f , and iiell. 

b) Simultaneously integrate forward in time to t^ the 

2n non-linear differential equations, Z = F(Z,t), and 

* 

the 2n linear perturbation equations, 6Z = 6Z. 

c) At the terminal time, evaluate h, llh|l, and lld|| to 
select the scalar correction factor n. 

d) Solve for 5X q and dt^. to update previous initial 
estimates. 

e) Examine II hi I and number of iterations to determine 
whether to continue iterating or stop. If the process 
is continued, return to step (a). 
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The following chapter will examine some alternate approaches 
to the solution of a TPBV problem which does not make use of the 
perturbation on h initially, Equation (4-4), but first considers some f 

function of h as a quantity to be minimized. Minimization of this \ 

function then results in a solution to the TPBV problem by switching 
to IMPF at some predetermined small value of h. This is done to take 
full advantage of the IMPF convergent characteristics which are very 
rapid whenever the initial guesses fall within the envelope of 
convergence. This technique may be classified as a direct-indirect 
approach; direct in the sense that some function of h is extremized 
directly, and indirect in the sense that IMPF does not extremize h 
to arrive at a solution. 






CHAPTER V 

MODIFIED METHOD OF PERTURBATION FUNCTIONS 
5.1 Introduction 

In the previous chapter, an improved method of perturbation 
functions (IMPF) was used to solve a TPBV problem arising out of the 
calculus of variations. This method essentially used the linear per- 
turbations of Z = F(Z,t), Equation (3-23), and the dissatisfaction 
vector h, Equation (3-24), to arrive at an initial correction vector 
, Equation (4-8b), with which to correct the previous initial X 

and terminal time estimates. The algorithm then makes these correc- 
tions and the process is continued as described in the last section of 
the previous chapter. Although this method is one of the most popular 
for solving TPBV problems, it is by no means independent of the initial 
estimates of co-states and terminal time. For this reason, three new 
strategies under the general name of modified method of perturbation 
functions (MMPF) are explored. These new strategies make use of 
function minimization methods which have shown good results in the 
literature and are used with the MPF method to form a direct-indirect 
approach to the solution of a TPBV problem. Essentially they are used 
to define an initial correction vector replacing Equation (4-8b). 

The name MMPF is derived from the fact that some of the features of 
MPF are used in MMPF. Most of the methods surveyed throughout the 
literature may be classified as either gradient methods, pattern 
search methods, or random search methods (an example is found in [44]). 
Of the three classifications, gradient methods generally use the most 
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intelligence and random search the least. The three particular 
methods chosen are (1) the Davidon Method [28] which Is a proven 
gradient technique; (2) the method of Hooke and Jeeves [26], a 
pattern search concept; and (3) an accelerated random search 
strategy [29]. The above three methods are identified as MMPF-D, 
MMPF-HJ, and MMPF-ARS, where the letter or letters after the hyphen 
in each case refers to Davidon, Hooke and Jeeves, and accelerated random 
search, respectively. For the sake of completeness, these individual 
methods will be briefly explained in the following three sections. 

In all cases of the modified method of perturbation 
functions, three changes had to be made to the method of perturbation 
functions concept described in Chapter IV to implement the direct- 
indirect concept of solution. The first change was the selection 
of a stopping condition on which to end the integration process 
since, with function minimization processes, there is no way of 
obtaining a terminal time estimate. This was accomplished by 
defining the terminal time as the time which satisfied one of the 
selected terminal boundary conditions in h; preferably, one which 
is monotonic decreasing. Therefore, the satisfaction of 

a = 0 (5-1) 

is used as the stopping condition. Next came the Silection of the 
function evaluated at the terminal time which is to be minimized. 

Since the TPBV problem is solved whenever the norm ^h reaches some 
small specified value, the logical choice for a function is some 
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U(h)>0. Thus, the function chosen to be minimized was 

U = h T h ^ 0 (5-2) 

It should be pointed out that the gradient information of 
U which is required in the MMPF-D would already be available from dh, 
Equation (4-4). Finally, since the set of the variables Z which 
minimize U are at the terminal time, they must be expressed as correc- 
tions to the initial Z's. This is accomplished by considering the 
solution of 6 X q from 6Z^ = 1T 5 ^ l5X 0 > Equation (4-7), that minimizes 
the sum of squares of the residuals. This solution is given by 


"o = <"« ’m*' 1 6Z f (5 ‘ 3) 

where 6Z^ is defined by the difference between Z at tp and the desired 
Z which results from the minimization of U. 

By making use of Equations (5-1) through (5-3), the MMPF 
!. |C ing one of the three previously mentioned function optimization 
methods would proceed until llhll reached some small predetermined value. 
At this time, a switch would be made to the IMPF strategy described in 
Chapter IV. A general computational outline of the direct- indirect 
method for solution using any of the three MMPF methods with the IMPF 
method as a final strategy is shown in Figure 4-1. 

5.2 Modified Method of Perturbation Functions Using the Davidon Method 
The method of Davidon as described by [28] involves the use 
of the first and second partial derivative information on the function 
being extremized and is derived assuming a quadratic function. 

Generally speaking, this method uses the idea of generating the 
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direction of search by multipling the gradient by a matrix that is in f j 

some manner related to the inverse of the second partial derivative 
matrix. This method appears particularly attractive to solve TPBV 
problems when used in conjunction with IMPF. Methods using gradient 

\ 

information (first partial derivatives) become extremely slow as the 

optimum is approached. Therefore, it is necessary to switch to a 

method which converges rapidly as the optimum is approached. As 

pointed out in Chapter IV, the IMPF procedure converges rjtpidly in 

the vicinity of the optimum, and a natural transition would be to ! 

this method. 

i'avidon's method described here briefly for a general 
function starts by examining a quadratic objective function, for 
example, 

U (Z) » a + b T Z + \ Z T cZ (5-4) 

where the .. ^ian of U(Z) (c) is a non-singular symmetric positive 
definite matrix. Tnus, U\Z) has a local minimum which, in fact, is 
also a global minimum. The gradient of (5-4) at any point Z is 

lf(Z) = b + cZ (5-5a) 

from which 

Z = c" 1 (iT(Z)-b) (5-5b) 

At the optimum of a differentiable function, all first derivatives 
are zero, and Equation (5-5b) may be written as 

Z* = -c -1 b (5-6) 
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The above equation and (5-5b) may be subtracted to obtain the correction 
A to be made to Z in order to have the optimum U, 

Z* - Z = A = -c _1 if (5-7) 

According to Davi don's method, if this inverse were available 
and the U under consideration were quadratic, the optimum could be 
found in one step. In most practical applications, this is not the 
case, and Equation (5-7), as such, is not applicable directly to most 
functions. In Davidon's method. Equation (5-7) is used only to define 
a direction of search in the space of the variables Z. And thus, a 
new correction vector, a, must be defined so that 

o = Sa = -S c' 1 if (5-8) 

where the step size, S, is found by obtaining the optimum along a line 
of search using quadratic interpolation. The evaluation of c"^ , as 
outlined in [28], is performed indirectly by using the change in 
gradients from iteration to iteration, the last c”^ , and the 
correction vector o. 

An iterative scheme for a minimum under MMPF-D may be 
described as follows: 

a) Begin with the identity matrix, I, as the initial guess 
for c," 1 , and assume some X Q . 

b) Integrate forward both the non-linear and the perturba- 
tion differential equations, as described in Section 4.4, 
observing Equations (5-1) through (5-3). 


i 


i 



1 


t 


,1 

.f 

E 

t 

42 

c) Evaluate the gradient of U, U . 

d) Compute a direction using Equation (5-7). 

I e) Evaluate the step size, S, by using a quadratic 

interpolation scheme or some other suitable scheme. 

f) Compute the correction vector, a, as 

g) Evaluate the desired terminal variables as a Z k + 
and compute 6\ using Equation (5-3). 

h) Compute the Qj^ for the next Iteration, which must be 
be positive definite so that the new direction will 
give local improvement of U(Z). 

i) If ilhll < 1/2 switch to IMPF; otherwise, start 

over at (b) with new estimates for A . 

5.3 Modified Method of Perturbation Functions Using the Method of 
Hooke and Jeeves 

The fiiethod of Hooke and Jeeves [26] is a direct pattern 
search technique requiring no derivative information but involving 
a sequential examination of trial solutions. This method is comprised 
of a series of exploratory moves about each variable In the function 
under consideration. These exploratory moves are performed on one 
variable at a time insisting that the value of the function be less 
than, or equal to, its previous value if a minimum is being sought. 

At the end of the total exploration in the variable space, it is 
assumed that a pattern may be established. Thus, a new exploration 
point is set by simultaneously stepping off a distance equivale •■j 
that of the initial start point and that at which the search stopped. 
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This technique is continued until the function reaches a desired 
value. The technique has strategies which handle situation where 
no improvement is detected whenever the establishment > - • rn 

is being sought. 

For the sake of completeness, an example of an applicaticn 
of this method in two-variable space (x,y) as pr posed by Hooke and 
Jeeves is presented. It is assumed that the search is startru at ‘he 
x and y corresponding to b-j in Figure 5-1, where b represents a 
pattern resulting from an exploratory search at E, and the subscripts 
indicate the trial. The resulting exploration about poin* ’ in 
Figure 5-1 becomes b£, and, by using b-| and and reasoning that a 
similar exploration about b^ would give similar results, a new 
exploration point, E 2 , is selected. Its location in space is given 
by 2 (b 2 - h-| ) ; using this formula, the pattern will continue to 
grow. The search continues until a point where improvement in the 
desired function is not realized (point E^ in Figure 5-1). At this 
point, a new strategy is needed to continue. The point E 4 is 
discarded, and the last good pattern point now becomes the new 
exploration point. Figure 5-2 shows that the new exolcration point 
has required a decrease in step to reach a new base bg. How using 
bg and bg, it is possible to start as before. The strategy continues 
until a solution is reached. 
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Figure 5-1. Schematic of Hooke and Jeeves Pattern 
Search Showing Improvement . 
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Figure 5-2. Schematic of Hooke and Jeeves Pattern Search 
Showing Pattern Destruction . 
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Extending the ideas presented by the example, an iterative 
scheme using IWPF-HJ for a minimum solution may be described as 
fol 1 ows : 

a) Integrate forward both the non-linear and the perturba- 
tion differential equations, as described in Section 4.4, 
observing Equations (5-1) and (5-2). 

b) Using the Hooke and Jeeves algorithm outlined in [26], 
sequentially compute each desired Z resulting from 
either an exploratory step or a pattern step. 

c) Compute SX Q according to Equation (5-3). 

d) If ||hil < 1/2 ilhil k _i switch to IMPF; otherwise, start 
over at (a) using the new estimates for x Q . 

It should be pointed out that using the method of Hooke and 
Jeeves requires at least 2n iterations before a pattern point is 
established. This is due to the fact that, during the exploratory 
phases, the effects of changes on each variable must be sequentially 
examined. 

5.4 Modified Method of Perturbation Functions Using an Accelerated 
Random Search MethotT " 

The accelerated random search method [29] is one which 
combines the ideas of both pattern search and random search with 
recalculation. This method has shown promise, particularly on 
multimodal hypersurfaces. The accelerated random search begins by 
making a random search in the vicinity of the initial estimate based 
on a normal distribution of probability. It follows with a step in 
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the opposite direction if the function is increased (failure), for 
example, in the case of search for a minimum. If this change In 
direction still produces no improvement, a new random step is 
generated (new set of . As soon as the random search phase 
discovers a direction for improvement (success), a pattern search 
is begun doubling the step after each success. However, if the 
pattern search produces a failure, a retreat is made to tne last 
improved step, end a random search is initiated. For practical 
reasons, the accelerated random search method confines its random 
search phase to an experimental region defined for each problem. 

The normal distribution function which governs the 
magnitude of the increments, 6Z, during the search phase is 



(5-9) 


where Y represents 2n random numbers uniformly distributed between 
0 and 1 and a is their standard deviations about a ^ean of zero. 

In addition, the sign of eZ^ has equal probability of being plus or 
minus. During the pattern phase, the increments assume an arithmetic 
progression form as 


5Z k+1 = 2 «Z k (5-10) 

and in the event of a failure during the random search phase, a 
reversal is executed by taking a step in the opposite direction as 

sZ w'- 24Z k < 5 - n ) 
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An iterative scheme using the HMPF-ARS now may be described 


as follows: 


a) Integrate forward both the non-linear and the 
perturbation differential equations, as described in 
Section 4.4, observing Equations (5-1) and (5-2). 

b) Using Equation (5-9), compute the 2n desired terminal 
corrections where, for each correction, a different Y 
is obtained. 

c) Transform the terminal corrections to initial 
corrections by using Equation (5-3) and proceed as in 
(a) with new initial x estimates. 

d) If U^ + i is an improvement over U k by either a random 
forward step (Equation 5-9), a reversal step 
(Equation 5-11), or a new random search, compute 

the next terminal correction by using Equation (5-10' 
and proceed as in (c). Continue the pattern search 
as long as there is improvement in U. Whenever 


llhli < 1/2 ||h|| k=1 , switch to IMPF. 

e) If there is no improvement in U during a pattern search, 
retreat to a random search mode, step (c), and start 
again. 

Both this method and the Davidon method make all the 
corrections simultaneously instead of sequentially as in the method 
of Hooke and Jeeves. 
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CHAPTER VI 

DISCUSSION AND ANALYSIS OF RESULTS 
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6.1 Introduction 

The numerical aspects involved in the solution of the TPBV 
problem will be discussed. In the main, this will include those which 
are common to all four methods including the boundary condi tic.is, 
numerical integration scheme used to integrate both the non-linear 
and linear differential equations, and finally the features which 
apply to each particular method. In this chapter, results are also 
shown which compare the IMPF and the three MMPF strategies outlined 
in Chapters 4 and 5, respectively, along with an additional strategy 
similar to that of [23]. Results are presented for the entry problem 
posed in Chapter II and transformed from an optimization problem into 
a TPBV problem as shown in Chapter III. 

It should be pointed out that all quantities appearing in 
this chapter, except the co-states, are dimensional only for the sake 
of awareness. Internally, all computations involved non-dimensional i zed 
quantities as stated earlier. 

6.2 Initial and Terminal Boundary Conditions 

The specification of the initial boundary conditions proceed 
as specified in Chapter III. First by selectng an initial state at 
t * t * o, an L (x , t ) vector may be stated as 
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L = 


= 0 


Z } (t ) - 25,975 ft/sec 
Z 2 (t Q ) - 400,000 ft 
Z 3 (t Q ) - 1.5 deg 
Z 4 (t 0 ) - 90 deg 
Z 5 (t Q ) - 90 deg 

[h (t o f - 0 deg j 

using the form of Equation (3.23). Likewise, a specified terminal 


(6-1) 


state vector may be selected at t = t^ as 


Z 1 (t^) - 25,306 ft/sec 
Z 2 (t f ) - 252,784 ft 

M a ^3 " 1.0568 deg _ q (6-2) 

Z 5 (t f ) - 89.9676 deg , 

Z 6 (t f ) - 18.1897 deg 


Equation (6-1) and the fact that the initial state and time are 
specified result in no initial transversal ity conditions from 
Equation (3-12). Thus, at t = t = o, the initial boundary condition 
becomes 


g = L = 0 


(6-3) 


On the other hand, since Equation (6-2) represents only five boundary 
conditions and seven are needed, the remaining two must come from the 
terminal transversal ity conditions, Equations (3-13) and (3-14), For 
the example problem, <j> = 0 and the terminal time t f and Z^ (t f ) have 
not been specified. Therefore, the remaining two boundary conditions 
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become H(t f ) = 0 and x 4 = ) a °* Thus, the terminal boundary 

condition h becomes 

Pz 1 (t f ) - 25,306 ft/sac"| 

Z 2 (t f ) - 252,784 ft/sec 
Z 3 (t f ) - 1.0568 deg 

h a Z 5 (t f ) - 89.9676 deg = 0 (6-4) 

Z 6 (t f ) - 18.1897 deg 

Z 10 ( V 

H (t f ) 

and as pointed out earlier it is the norm of h, jlhil , that is used to 
determine convergence to a solution. Furthermore, since time does not 
appear explicitly in the generalized Haiultonian 


H (t) = 0 (6-5) 

for t Q < t £ t^ where H is defined using Z variables as 

H _ r 7 3.15 -0.5 .1 A lr 2 r 2v 0 * 5 - .2 Jf 6 7 F ,, ,, 

H " C 1 Z 1 p + 7 Z (C L + C D } p Z 1 + .f Z i+6 F i (6 " 6) 

W 1 “ I 


For this particular problem, the control, u, is eliminated 
from the Z equations by using the optimality condition, H u = 0, and the 
Legendre-Clebsch condition for a minimum, H uu > 0, as outlined in 
Section 3.2. Use of the optimality condition yields 


x 3 sin u 


cos u 
a 4 cos r 


(6-7) 


which may be used to solve for the -functions sin u and cos u by 
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squaring (6-7) and using the trigonometric identity sin u+cos u = 1 
to give 

*4 

sin u = (6-8a) 

V ~2 2 ? — 

A 4 + A3 cos r 


cos r 

cos u = — _ ■ 
+ Vx 4 + cos 2 r 


(6-8b) 


The Legendre-Clebsch condition is then used to choose the proper sign 
for Equation (6-8). This condition requires that 


H uu = *3 cos u "*4 MtT 


> 0 


(6-9) 


and using Equation (6-8) to substitute for sin u and cos u in 
Equation (6-9) shows that for -iT/2<r<ir/2 the negative sign is chosen 
for Equation (6-8a) and the positive sign for Equation (6-8b) or 


sin u = - (6-10a) 

V r^ 2 7 — 

x 4 + x^ cos r 

x, cos r 

cos u = - (6-1 Ob) 

It should be pointed out that, as a result of the terminal trans- 
versality conditions (6-4), x^(t^) = 0 and therefore u(t^) = 180 degrees 
on the optimal trajectory. 
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6.3 Numerical Integration Method 

The numerical integration method used [45,46] to integrate 
the 2n non-linear differential equations, Z =■ F{Z,t), and the 2n linear 

• aC 

perturbation equations, 67. = 6Z, is a variable step, fourth-order 

• aF 

Runga-Kutta scheme. The Z equations of motion and the ^ are 

functionally presented in Appendix A. Since the entry problem chosen 

-as an example is not for fixed terminal time, the generalized 

Hamiltonian becomes time-independent and, therefore, constant. This 

fact is used as an aid in selecting relative error bounds for the 

variable time step integration process. The task of selecting error 

bounds, particularly upper bounds, is a formidable one and often can 

be selected only according to experience and the type of problem at 

hand. The fact that the generalized Hamiltonian is a constant, and 

zero on the converged solution, allows observation of the effects of 

the error bound on the numerical accuracy and thus makes possible a 

logical choice. Table 6-1 shows the effects of the error bound for 

six cases on the variation of the generalized Hamiltonian during an 

iteration, as well as the value of the norm at the end of the 

iteration and the required UNIVAC 1108 time. The method of solution 

used was the IMPF described in Chapter IV; all six cases used the 

same identical boundary conditions and final time estimate. It is 

logical, therefore, to assume that this table gives a representative 

relative assessment of the effects of error bound on the accuracy 

of solution. Based on this information, a relative error bound of 
-3 -5 

5.0 x 10 - 5.0 x 10 was chosen based on a convergence criterion 
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The Influence of Relative Error Bound Magnitude with an Error Bound 
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of h * 1x10 compared to an expected Hamiltonian amplitude of about 
the same order of magnitude. A further consideration in the choice of 
the error bound was the exponential increase nature of the computer 
time with decrease in error bound. 

6.4 Solution of the Entry Problem by IMPF 

Solution of the entry problem using the terminal conditions 
of Section 6.2 and the integration method and constants stated in 
Section 6.3 was initiated by assuming a terminal time, t^, of 260 
seconds with which to begin the iteration process outlined in 
Chapter IV. The necessary constants needed for the correction process 
as described in Section 4.3 are as follows for jidil < lie l! «i/7i 


f0.4 if hhli > 1 x 10 
1 0. 4 if ilh|l < 1 x 10‘ 
1 0. 8 if ilhil < 1 x 10' 
0.0 if ilhil < 1 x 10“ 


where h is assumed to be evaluated using non-dimensional i zed 
variables. The choice of |te|l= v7 comes from allowing an absolute 
range of unity for all the states and time since they are all non- 
dimensionalized and result in magnitudes close to one. The vector 

* 3 h • 

h and the elements that are necessary to compute the corrections 
according to Equation (4-8a) are presented in Appendix A. 

6.5 Solution of Entry Problem by MMPF 

The stopping condition, ft, chosen to end the integration 
process for each iteration in the MMPF solution was 

ft = Z 2 (t f ) - 252,784 ft = h 2 = 0 (6- 
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This variable was monotonic decreasing for all the c?ses studied and, 
therefore, served as an excellent stopping condition. The function 
chosen to be extremized resulted in a highly non-linear function 
involving all the stites and co-states as well as the state differential 
equations. This function may be written as 

U * h^h = h^ + h^ + h^ + h^ + hg + hg + hy (6-13) 

where h is defined by Equation (6-4). The value of llhll at which the 
transition from MMPF to IMPf was executed was set at iihll = -* 2 “ I k=l . 

This particular value was chosen under the assumption that llhll = k=l 

is sufficiently close to the minimum so that the IMPF can attempt to 
obtain a solution. Thus, at the above value of ilhll, the IMPF is 
started with not only the co-state values o ; the particular iteration 
but also with the corresponding final time as the terminal time estimate. 
It can be expected that whenever the transition to IMPF is made, the 
final time estimate will Le different than t^=260 seconds, which was 
assumed to start the IMPF solutions of Section 6.4. 

6.6 Results of Solution of ;he Entry Problem by IMPF and MMPF 

In order to adequately explore the potential of the IMPF and 
the MMPF to solve the TPBV problem, a set of six initial co-state (x) 
vectors wer.: selected as shown in Table 6-2. The selection of the 
first four vectors made use of the information from the terminal 
boundary conditions which required the control u to be positive; thus, 
the first three co-states were chosen positive and the fourth, negative. 
Since it is apparent from Table 6-2 that the choice of the first vector 
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is close to the optimal, the choice of the remaining five co-states was 
made to measure the effect on convergence of starting further away 
from the optimum. To further evaluate the improvement of these two 
methods, an additional MPF strategy based solely on the magnitude of 
the norm and similar to that of Lewallen [23] was used. The data in 
Table 6-2 show results for the following examples: Cases 1 through 6 

used the IMPF with its correction factor based on norm and step size; 
Cases 7 through 10 used the MPF with a correction factor based only on 
norm and similar to that of Equation (6-11 ) with the exception that 
the first factor is 0.2 instead of 0.4; and finally, Cases 11 through 
13 used the MMPF where a, b, and c stand for Davidon, Hooke-Jeeves, 
and accelerated random search, respectively, with a transition to the 
IMPF of Cases 1 through 6 as explained in Section 6.5. 

The IMPF strategy using a correction factor dependent on the 
norm and the step size was used to solve the TPBV problem starting with 
each of the six co-state vectors of Table 5-2 as initial estimates and 
guessing an initial terminal time t^ = 260 seconds. This strategy was 
successful in obtaining a solution for each of the six candidate 
co-state vectors. The results depicted in Table 6-3 show essentially 
that the vector closest to the optimal solution required the least 
number of iterations. Cases 2 through 4, which used somewhat more 
difficult vectors, also converged in essentially the same number of 
iterations. Cases 5 and 6 are by far the most difficult but yet were 
driven to convergence by the strategy. In summary, the solution of 
the TPBV problem was very insensitive to the first four cases. However, 
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as the progression was made to order of magnitude differences in the 
initial co-state (Case 6), the sensitivity increased, but the most 
important fact was that a solution was still obtained. This point is 
further emphasized when the MPF with a correction factor dependent on 
the norm only is discussed next. 

Results using a correction factor philosophy similar to [23] 
and dependent on the norm only are presented in Table 6-4. As shown 
on Table 6-2, this scheme used only the first four vectors in Cases 7 
through 10, respectively. The only solution possible out of the four 
cases was for Case 7 which required 50 percent more iterations than 
Case 1 for the same co-state vector and terminal time estimate. The 
remaining three cases diverged after four iterations even though the 
first two iterations gave a decrease from the first norm. This set 
of cases points out the power of a correction factor such as the one 
used in Cases 1 through 6. 

Finally, results using the MMPF and starting with the first, 
fourth, and fifth co-state vectors are depicted in Table 6-5. The 
logic used to control the total number of iterations for these cases 
was as follows. A maximum of 50 iterations was allowed within which 
to halve the initial norm; once the transition was made to IMPF, a 
maximum of 30 iterations was allowed if the total number of iterations 
exceeded the companion case of Table 6-3. However, if the total 
number was less, iterations were continued up to an equivalent total 
as tne companion case of Table 6-3. Table 6-5 exhibits two cases, 12a 
and 13c, which converged. Although Case 12a converged, it did so in 
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about one-third more iterations than Case 4. Case 13c, on the other 
hand, required approximately 26 percent fewer iterations than its 
companion case, Case 5. The remaining cases shown in Table 6-5 did 
not converge, nor did they diverge In the set number of Iterations. 

It is felt that the cases which were evaluated will converge if allowed 
sufficient iterations but there was no reason for doing this. Based 
on the cases examined in Table 6-5, the MMPF using the accelerated 
random method seems the most powerful and it appears to get relatively 
better as it starts with co-states further from their optimum, as 
evidenced by Cases 11c, 12c, and 13c. The method of Hooke-Oeeves does 
not appear promising in comparison. Davidon's method converged on 
Case 12a but did not meet the basic criteria in 11a and 13a. This is 
not surprising because gradient methods are incapable of negotiating 
multimodal hypersurfaces and Equation ( 6 - 13 ), as evidenced by examina- 
tion of Cases Ua and 13a, appears to be multimodal. 

Figures 6-1 and 6-2 depict the optimum states, control 
variable, and the pay-off function variables Q and G as a function 
of the optimal time. It should be pointed out that, although the 
theory only guarantees a local optimum, all of the cases which 
converged do so to the identical solution. Thus, the statement may 
be made that the solution obtained and shown in Figure 6-1 is the 
optimum solution at least within the space defined by the six co-state 
vectors of Table 6-2. 
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Figure 6-1. Optimal States for Solution of Entry Problem 
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CHAPTER VII 

SUMMARY AND CONCLUSIONS 

Two general iterative methods employing the basic method of 
perturbation functions to solve TPBV problems have been described. 

The first method ( IMPF ) uses a unique selection scheme to arrive at 
a correction factor which in essence sets the desired changes to be 
applied to the previous guessed co-states. The second method (MMPF) 
solves the problem in two phases. Phase 1, the start, uses function 
minimization techniques until the norm is subsequently reduced to a 
level equal to, or less than half of, that which it began with. At 
this point Phase 2 which is identical to the first method, takes over. 
To examine the capabilities of these two methods to obtain a solution, 
a complex entry optimization problem which was reduced to a TPBV 
problem through calculus of variations techniques was used. The two 
methods were then exercised using the above problem starting with 
several sets of different initial co-states. In addition, a method 
similar to that used in [23] on an Earth-Mars transfer problem was used 
to serve as an additional basis of comparison. 

Based on this study, both of the techniques used must be 
judged superior to the one which is similar to that of [23]. Of the two 
techniques used, the IMPF with the correction factor based on norm and 
step size is by far superior as evidenced by its ability to obtain a 
solution for all the candidate cases. However, the MMPF using the 
accelerated random search strategy seems to become more attractive as 
the initial co-states get further away from their optimum values. Due 
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to the extreme complexity of the problem considered, it became 
computationally impractical to study cases further away from the 
optimum than those attempted. 

Though convergent properties of iteration schemes are problem 
dependent, certain observations may be made in comparing the results of 
this study to those in the literature. Tapley and Lewallen [41] 
presented an Earth-Mars transfer problem which exhibited convergence 
from regions much closer to the optimum than those presented here. A 
similar statement may be made of studies on reentry problems such as 
the two-dimensional problem studied by Lastman and Tapley [19]. A 
recent study by Tapley and Williamson [32] on a three-degree-of- freedom 
reentry problem using MPF indicated 15 iterations were required for an 
optimal solution. However, based on the close proximity of their 
initial co-state estimates, indications are that the IMPF technique 
might improve on the iterations required for an optimal solution. In 
conclusion, it may be said that two methods have been formulated which 
have decreased the contingency of a solution to a complex entry problem 
on the initial co-state estimates. 
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APPENDIX A 
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TH E STATE, CO-STATE, AND MPF EQUATIONS 

The state (x) and co-state (x) differential equations of 
motion for the entry problem form a set of 2n = 12 equations which are 
ncn-linear and ordinary in nature. Of this set of 2n equations, n 
are state equations and n are co-state equations. Both sets of equations 
may be expressed as functions of the variational Hamiltonian. The 
variational Hamiltonian (Equation (6-6)) can be written in terms of , 
and u variables as 


H = C,x 


1 A 1 


3.15 “0.5 . 
P + 
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4 (c^ + 
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r 2J /2 “ 2 y 

Cp. ) p x-i + 2-i 

u 1 i=l 


x i f i 


(A-l ) 


where the n set of states x correspond to Equations (2-3a) through (2-3f) 
respectively; p(x^) is of the form explained in Section 2.2; and 
F(x,u,t) represents the n differential equations of motion (2-3a) 
through (2-3f). Taking the partial derivative of H with respect to x 
yields n differential equations of state 

x = H x = F(x,u,t) (A-2) 

The n equations of co-state (x) are then obtained by taking the partial 
derivative of H with respect to x 

X = H x = x(x,X,u,t) (A-3) 

Equations (A-2) and (A-3) may be combined to form a 2n set of differential 
equations 

Z = F(Z,u,t) (A-4) 

where the variable Z is made up of the x and x variables. 
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The perturbation equations, result in a 2n X 2n matrix 
which may also conveniently be expressed as a function of H 


3F . V 1 ) 


D 3 (t) 


D 2 (t) 




v*’ ■ H xx - h xu h ;1 h . 
D 2<‘> ' H xu H ux 


°3<‘> ■ - H xx * H xu H u x < A - 8 > 

For example, with the entry problem defined in this thesis, the top row 
of the partition is as follows, with u being a scalar since there is 
but one control variable: 
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Referring to Equation (6-4), the vector equation for the time 
derivatives of the terminal dissatisfaction h is 
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where H(t f ) = 0 since H(t) is a constant, as pointed out in Section 6.2, 
and Z(t f ) are the 2n differential equations of state and co-state. The 
results in a 7 x 12 matrix which has zero for many of its elements. 
Only the non-zero elements will be specifically identified 
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The last row of the matrix is which is 
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